home *** CD-ROM | disk | FTP | other *** search
/ Technotools / Technotools (Chestnut CD-ROM)(1993).ISO / lang_oth / linpkdrv / sgt.for < prev    next >
Text File  |  1985-01-12  |  9KB  |  321 lines
  1. C     MAIN PROGRAM
  2.       INTEGER LUNIT
  3. C     ALLOW 5000 UNDERFLOWS.
  4. C     CALL TRAPS(0,0,5001,0,0)
  5. C
  6. C     OUTPUT UNIT NUMBER
  7. C
  8.       LUNIT = 6
  9. C
  10.       CALL SGTTS(LUNIT)
  11.       STOP
  12.       END
  13.       SUBROUTINE SGTTS(LUNIT)
  14. C     LUNIT IS THE OUTPUT UNIT NUMBER
  15. C
  16. C     TESTS
  17. C        SGTSL,SPTSL
  18. C
  19. C     LINPACK. THIS VERSION DATED 08/14/78 .
  20. C     JACK DONGARRA, ARGONNE NATIONAL LABORATORY.
  21. C
  22. C     SUBROUTINES AND FUNCTIONS
  23. C
  24. C     LINPACK SGTSL,CTPSL
  25. C     EXTERNAL SGTXX
  26. C     BLAS SASUM
  27. C     FORTRAN ABS,AIMAG,AMAX1,FLOAT,REAL
  28. C
  29. C     INTERNAL VARIABLES
  30. C
  31.       INTEGER LUNIT
  32.       REAL B(20),BSAVE(20),D(20),EYE,C(20),E(20),X(20)
  33.       REAL ANORM,EN,ENORM,EPS,Q(2),RNORM,SASUM,XNORM
  34.       REAL SMACH
  35.       INTEGER I,INFO,IPT,PD,KASE,KFAIL(2),KSING,N,NM1,NPRINT,POSDEF
  36.       EYE = 0.0E0
  37. C     REAL EYE = IMAGINARY UNIT, REAL EYE = ZERO
  38. C
  39. C
  40. C     WRITE MATRIX AND SOLUTIONS IF  N .LE. NPRINT
  41. C
  42.       NPRINT = 3
  43. C
  44.       WRITE (LUNIT,230)
  45.       DO 10 I = 1, 2
  46.          KFAIL(I) = 0
  47.    10 CONTINUE
  48.       KSING = 0
  49. C
  50. C     COMPUTE MACHINE EPSILON
  51. C
  52.       EPS = SMACH(1)
  53.       WRITE (LUNIT,240) EPS
  54.       WRITE (LUNIT,220)
  55. C
  56. C     START MAIN LOOP
  57. C
  58.       KASE = 1
  59.    20 CONTINUE
  60. C
  61. C        GENERATE TEST MATRIX
  62. C
  63.          CALL SGTXX(C,D,E,N,KASE,POSDEF)
  64. C
  65. C        N = 0 SIGNALS NO MORE TEST MATRICES
  66. C
  67. C     ...EXIT
  68.          IF (N .LE. 0) GO TO 210
  69.          INFO = 0
  70.          PD = 1
  71.          IF (POSDEF .EQ. 1) PD = 2
  72.          DO 200 IPT = 1, PD
  73.             WRITE (LUNIT,250) KASE
  74.             WRITE (LUNIT,260) N
  75.             IF (N .GT. 1) GO TO 30
  76.                ANORM = ABS(D(1))
  77.                WRITE (LUNIT,450) D(1)
  78.                X(1) = 1.0E0
  79.                B(1) = D(1)
  80.                BSAVE(1) = B(1)
  81.             GO TO 110
  82.    30       CONTINUE
  83.                NM1 = N - 1
  84.                ANORM = ABS(D(1)) + ABS(C(2))
  85.                IF (N .LE. 2) GO TO 50
  86.                   DO 40 I = 2, NM1
  87.                      ANORM = AMAX1(ANORM,
  88.      *                             ABS(C(I+1))+ABS(D(I))+ABS(E(I-1)))
  89.    40             CONTINUE
  90.    50          CONTINUE
  91.                ANORM = AMAX1(ANORM,ABS(E(N-1))+ABS(D(N)))
  92.                WRITE (LUNIT,430) ANORM
  93. C
  94.                IF (N .GT. NPRINT) GO TO 60
  95.                   WRITE (LUNIT,220)
  96.                   WRITE (LUNIT,450) (C(I), I = 2, N)
  97.                   WRITE (LUNIT,220)
  98.                   WRITE (LUNIT,450) (D(I), I = 1, N)
  99.                   WRITE (LUNIT,220)
  100.                   WRITE (LUNIT,450) (E(I), I = 1, NM1)
  101.                   WRITE (LUNIT,220)
  102.    60          CONTINUE
  103. C
  104. C              GENERATE EXACT SOLUTION
  105. C
  106.                X(1) = 1.0E0
  107.                IF (N .GE. 2) X(2) = EYE
  108.                IF (N .LE. 2) GO TO 80
  109.                   DO 70 I = 3, N
  110.                      X(I) = -X(I-2)
  111.    70             CONTINUE
  112.    80          CONTINUE
  113. C
  114. C              SAVE MATRIX AND GENERATE R.H.S.
  115. C
  116.                B(1) = D(1)*X(1) + E(1)*X(2)
  117.                BSAVE(1) = B(1)
  118.                IF (N .LE. 2) GO TO 100
  119.                   DO 90 I = 2, NM1
  120.                      B(I) = C(I)*X(I-1) + D(I)*X(I) + E(I)*X(I+1)
  121.                      BSAVE(I) = B(I)
  122.    90             CONTINUE
  123.   100          CONTINUE
  124.                B(N) = C(N)*X(N-1) + D(N)*X(N)
  125.                BSAVE(N) = B(N)
  126.   110       CONTINUE
  127. C
  128. C           FACTOR AND SOLVE A GENERAL TRIDIAGONAL SYSTEM
  129. C
  130.             IF (IPT .EQ. 1) CALL SGTSL(N,C,D,E,B,INFO)
  131. C
  132. C           TEST FOR SINGULARITY
  133. C
  134.             IF (INFO .EQ. 0) GO TO 120
  135.                WRITE (LUNIT,270)
  136.             GO TO 190
  137.   120       CONTINUE
  138. C
  139. C              FACTOR AND SOLVE A POSITIVE DEFINITE SYSTEM
  140. C
  141.                IF (IPT .EQ. 2) CALL SPTSL(N,D,E,B)
  142.                IF (IPT .EQ. 1) WRITE (LUNIT,280)
  143.                IF (IPT .EQ. 2) WRITE (LUNIT,290)
  144.                IF (N .GT. NPRINT) GO TO 130
  145.                   WRITE (LUNIT,300)
  146.                   WRITE (LUNIT,460) (B(I), I = 1, N)
  147.                   WRITE (LUNIT,310)
  148.                   WRITE (LUNIT,460) (BSAVE(I), I = 1, N)
  149.   130          CONTINUE
  150. C
  151. C              COMPUTE ERRORS AND RESIDUALS
  152. C                 E  =  X - X
  153. C                 R  =  B - A*X
  154. C
  155.                XNORM = SASUM(N,X,1)
  156.                CALL SGTXX(C,D,E,N,KASE,POSDEF)
  157.                IF (N .GT. 1) GO TO 140
  158.                   RNORM = ABS(D(1)*B(1)-BSAVE(1))
  159.                   ENORM = ABS(B(1)-X(1))
  160.                GO TO 170
  161.   140          CONTINUE
  162.                   ENORM = ABS(B(1)-X(1))
  163.                   RNORM = ABS(D(1)*B(1)+E(1)*B(2)-BSAVE(1))
  164.                   IF (N .LE. 2) GO TO 160
  165.                      DO 150 I = 2, NM1
  166.                         RNORM = RNORM
  167.      *                          + ABS(C(I)*B(I-1)+D(I)*B(I)+E(I)*B(I+1)
  168.      *                                -BSAVE(I))
  169.                         ENORM = ENORM + ABS(B(I)-X(I))
  170.   150                CONTINUE
  171.   160             CONTINUE
  172.                   RNORM = RNORM + ABS(C(N)*B(N-1)+D(N)*B(N)-BSAVE(N))
  173.                   ENORM = ENORM + ABS(B(N)-X(N))
  174.   170          CONTINUE
  175. C
  176.                WRITE (LUNIT,320) ENORM
  177.                WRITE (LUNIT,330) RNORM
  178. C
  179. C              COMPUTE TEST RATIOS
  180. C
  181.                EN = FLOAT(N)
  182.                Q(1) = RNORM/(EPS*ANORM*XNORM)
  183.                Q(2) = ENORM/(EPS*XNORM)
  184.                WRITE (LUNIT,220)
  185.                WRITE (LUNIT,340)
  186.                WRITE (LUNIT,220)
  187.                WRITE (LUNIT,400)
  188.                WRITE (LUNIT,410)
  189.                WRITE (LUNIT,420)
  190.                WRITE (LUNIT,220)
  191.                WRITE (LUNIT,440) (Q(I), I = 1, 2)
  192.                WRITE (LUNIT,220)
  193.                IF (N .EQ. 1) EN = 2.0E0
  194.                DO 180 I = 1, 2
  195.                   IF (Q(I) .GT. EN) KFAIL(I) = KFAIL(I) + 1
  196.   180          CONTINUE
  197.   190       CONTINUE
  198. C
  199.             WRITE (LUNIT,350)
  200.   200    CONTINUE
  201.          KASE = KASE + 1
  202.       GO TO 20
  203.   210 CONTINUE
  204. C
  205. C     FINISH MAIN LOOP
  206. C
  207. C     SUMMARY
  208. C
  209.       WRITE (LUNIT,360)
  210.       KASE = KASE - 1
  211.       WRITE (LUNIT,370) KASE
  212.       WRITE (LUNIT,380) KSING
  213.       WRITE (LUNIT,390) KFAIL
  214.       WRITE (LUNIT,470)
  215.       RETURN
  216. C
  217. C     ALL FORMATS
  218. C
  219.   220 FORMAT (1H )
  220.   230 FORMAT (29H1LINPACK TESTER, SGT**, SPT** /
  221.      *        30H THIS VERSION DATED 08/14/78 .)
  222.   240 FORMAT (18H MACHINE EPSILON =, 1PE13.5)
  223.   250 FORMAT (14H MATRIX NUMBER, I4)
  224.   260 FORMAT (4H N =, I4)
  225.   270 FORMAT (16H MAYBE SINGULAR.)
  226.   280 FORMAT (18H RESULTS FOR SGTSL)
  227.   290 FORMAT (18H RESULTS FOR SPTSL)
  228.   300 FORMAT ( / 4H X =)
  229.   310 FORMAT ( / 4H B =)
  230.   320 FORMAT (14H ERROR NORMS =, 1P2E13.5)
  231.   330 FORMAT (14H RESID NORMS =, 1P2E13.5)
  232.   340 FORMAT (26H TEST RATIOS.. E = MACHEPS)
  233.   350 FORMAT ( / 14H ************* /)
  234.   360 FORMAT (8H1SUMMARY)
  235.   370 FORMAT (18H NUMBER OF TESTS =, I4)
  236.   380 FORMAT (30H NUMBER OF SINGULAR MATRICES =, I4)
  237.   390 FORMAT (30H NUMBER OF SUSPICIOUS RATIOS =, 8I4)
  238.   400 FORMAT (20H    ERROR     RESID )
  239.   410 FORMAT (2(10H   -------))
  240.   420 FORMAT (20H     E*X      E*A*X )
  241.   430 FORMAT (14H NORM(A)     =, 1PE13.5)
  242.   440 FORMAT (8F10.4)
  243.   450 FORMAT (6G11.4)
  244.   460 FORMAT (2G14.6)
  245.   470 FORMAT ( / 12H END OF TEST)
  246.       END
  247.       SUBROUTINE SGTXX(C,D,E,N,KASE,POSDEF)
  248. C     FORTRAN FLOAT
  249. C
  250.       INTEGER N,KASE,POSDEF
  251.       REAL C(1),D(1),E(1)
  252. C
  253.       REAL EYE
  254.       INTEGER I
  255. C
  256.       EYE = 0.0E0
  257.       GO TO (10,20,30,30,30,50,50,70,70,90,110), KASE
  258. C
  259.    10 CONTINUE
  260.          N = 1
  261.          D(1) = 1.0E0
  262.          POSDEF = 1
  263.       GO TO 120
  264. C
  265.    20 CONTINUE
  266.          N = 2
  267.          D(1) = 4.0E0
  268.          D(2) = 4.0E0
  269.          C(2) = 2.0E0
  270.          E(1) = 2.0E0
  271.          POSDEF = 1
  272.       GO TO 120
  273. C
  274.    30 CONTINUE
  275.          N = (KASE - 2)*3
  276.          DO 40 I = 1, N
  277.             C(I) = 1.0E0/(FLOAT(2*I+2) + EYE)
  278.             D(I) = 1.0E0/(FLOAT(2*I+1) + EYE)
  279.             E(I) = 1.0E0/(FLOAT(2*I) + EYE)
  280.    40    CONTINUE
  281.          POSDEF = 0
  282.       GO TO 120
  283. C
  284.    50 CONTINUE
  285.          IF (KASE .EQ. 6) N = 10
  286.          IF (KASE .EQ. 7) N = 20
  287.          DO 60 I = 1, N
  288.             C(I) = 1.0E0
  289.             D(I) = 4.0E0
  290.             E(I) = 1.0E0
  291.    60    CONTINUE
  292.          POSDEF = 1
  293.       GO TO 120
  294. C
  295.    70 CONTINUE
  296.          IF (KASE .EQ. 8) N = 10
  297.          IF (KASE .EQ. 9) N = 20
  298.          DO 80 I = 1, N
  299.             C(I) = 1.0E0 + EYE
  300.             D(I) = 4.0E0 + EYE
  301.             E(I) = 1.0E0 - EYE
  302.    80    CONTINUE
  303.          POSDEF = 1
  304.       GO TO 120
  305. C
  306.    90 CONTINUE
  307.          N = 10
  308.          DO 100 I = 1, N
  309.             C(I) = 0.0E0
  310.             D(I) = 1.0E0
  311.             E(I) = 0.0E0
  312.   100    CONTINUE
  313.          POSDEF = 1
  314.       GO TO 120
  315. C
  316.   110 CONTINUE
  317.          N = 0
  318.   120 CONTINUE
  319.       RETURN
  320.       END
  321.